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Abstract 


A fundamental prediction of relativistic cosmologies is that, due to 
the expansion of space, observations of the distant cosmos should 
be time dilated and appear to run slower than events in the local 
universe. Whilst observations of cosmological supernovae unambigu- 
ously display the expected redshift-dependent time dilation, this has 
not been the case for other distant sources. Here we present the 
identification of cosmic time dilation in a sample of 190 quasars 
monitored for over two decades in multiple wavebands by assessing 
various hypotheses through Bayesian analysis. This detection counters 
previous claims that observed quasar variability lacked the expected 
redshift-dependent time dilation. Hence, as well as demonstrating the 
claim that the lack of the redshift dependence of quasar variabil- 
ity represents a significant challenge to the standard cosmological 
model, this analysis further indicates that the properties of quasars 
are consistent with them being truly cosmologically distant sources. 
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1 Introduction 


A fundamental consequence of the relativistic picture of expanding space is 
cosmological time dilation, where events in the distant universe appear to run 
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slowly compared to those in the local cosmos [1-3]. Whilst this time dilation 
has been unambiguously detected in the light curves exhibited by cosmo- 
logically distant supernovae [4-8], the appearance of time dilation in other 
cosmic sources is less conclusive. For example, whilst examinations of the light 
curves of gamma-ray bursts (GRBs) have generally shown consistency with the 
expected cosmological signature, uncertainties in the detailed emission mech- 
anism and expected light curve characteristics mean that this detection has 
not been definitive [e.g. 9-13]. Furthermore, the role of the more recently dis- 
covered fast radio bursts [FRBs: 14] as ’standard clocks’ is similarly limited 
by knowledge of the physical processes driving the output [15]. 

Quasars have been known to be variable sources since their discovery in the 
1960s [16], with emission arising from a relativistic accretion disk orbiting a 
supermassive black hole [17]. However, it has been claimed that the variability 
displayed by quasars over a broad range of redshifts does not show the expected 
cosmological time dilation [18-20]. This has led to the suggestion that quasar 
variability is not intrinsic, but is due to microlensing due to the presence of 
cosmologically distributed black holes [21, 22]. Others have stated that this 
points to more fundamental issues with our cosmological ideas [e.g. 23-25], 
with even the suggestion that quasars are not cosmologically distant and that 
their observed redshifts are due to mechanisms other than the expansion of 
space. 

In 2012, a study of the variability characteristics of a sample of thir- 
teen quasars observed behind the Magellanic Clouds as part of the MACHO 
microlensing program was suggestive of the expected (1 + z) time dilation 
dependence, where z is the quasar redshift [26]. However, with the small sam- 
ple and relatively short monitoring period, this result is inconclusive. Recently, 
a new sample of the variability properties of quasars was presented as part 
of the Dark Energy Survey and comprises 190 quasars, covering the redshift 
range z ~ 0.2 > 4.0 [27]. These are drawn from the sample of more than one 
hundred thousand spectroscopically identified quasars with absolute magni- 
tude M < —22 in the 300 deg? Sloan Digital Sky Survey of Stripe 82 (S82). 
Published as part of the SDSS DR7 quasar catalogue, the physical properties 
of these quasars are presented in [28]. This includes the bolometric luminosity 
which was determined through spectral fitting and correction from compos- 
ite spectral energy distributions [29]. These quasars were photometrically 
observed between 1998 and 2020, and so for more than two decades, through 
the combination of multiple epochs of exposures from SDSS, PanSTARRS-1, 
and the Dark Energy Survey, with additional follow-up monitoring with Blanco 
4m/DECam. 

The total dataset consists of roughly two hundred photometric observations 
of each quasar in multiple bands, although the cadence of these observations 
is very uneven over the observing period. To account for this when calculating 
characteristic timescales of the quasar variabilities, [27] adopted a Gaussian 
process regression [e.g. 30] to interpolate the photometric data and the associ- 
ated uncertainties between the observations; details are given in Appendix A2 
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of their paper. Each quasar light curve in each band is represented as a Damped 
Random Walk [DRW; 31, 32]; this is found to be an accurate description of 
quasar variability with only a mild dependence on the physical properties of 
the quasars [33]. Practically, this defines the covariance matrix of the Gaussian 
Process that describes the variability. With this, the Gaussian Process regres- 
sion software, Celerite [34], is used to determine the characteristic DRW time 
scale, as well as the 16” and 84*” percentiles of the distribution. Armed with 
these bolometric luminosities and variability time scales drawn from the DRW 
analysis, the goal of this paper is to search for the signature of cosmological 
time dilation of these distant sources. 


2 Results 


In the following analysis, we consider the redshift dependence of time dilation 
to be of the form (1 + z)”, where z is the redshift of the source. Clearly, for 
the expected cosmological dependence, n = 1, whilst n = 0 demonstrates no 
redshift dependence, representative of the claims from several previous studies 
of quasar samples [e.g 20]. To explore the various possibilities, several distinct 
hypotheses are explored. These are: 


e Ho: n is fixed at zero, representing no redshift dependence on the observed 
quasar timescales. 

e Hı: n is fixed at one, representing the expected redshift dependence of the 
cosmological time dilation. 

e Hə: n is treated as a free parameter. 

e Hs and Ha: n is fixed at —1 and 2 respectively. 


The final two cases represent extreme cases where additional influences, such 
as quasar evolution, may significantly influence the observed time variability 
of quasars. 

As outlined in the Methodology (Section 3), these differing hypotheses were 
compared through the calculation of the Bayesian evidence [35] for each situ- 
ation under consideration, with the results of these calculations presented in 
Table 1; in assessing the ratio of Bayesian evidences, a factor of 10-100 is con- 
sidered a strong favouring of one hypothesis over another, whereas greater than 
100 is decisive [36]. One immediate conclusion is that the favoured hypothesis 
is Hı, the case where n = 1, which represents the expected redshift depen- 
dence of the cosmological time dilation. This is significantly favoured over the 
alternative Ho, with an evidence ratio greater than 10°, which represents the 
situation where there is no redshift dependence on the observed timescales of 
cosmological variability. Furthermore, 1, is significantly favoured over the two 
extreme cases, H3 and H4. 

The posterior distribution for the redshift dependence for the time dila- 
tion, n, specifically H2, where this is treated as a free parameter, is presented 
in Figure 1. Reflecting the previous analysis, the is clearly offset from zero, 
indicative of a redshift dependence of the observed timescale of variability 
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over the quasar sample. This posterior distribution, which may be summarised 
using n = 1.28+0-29, is consistent with the expected cosmological dependence 
with n = 1, and the presented analysis significantly favours the presence of 


cosmological time dilation of the observed quasar variability. 


3 Methodology 


Probing the fundamental nature of our universe often calls upon standard 
rulers or candles to allow us to determine the influence of expansion on observ- 
able quantities. In hunting for cosmological time dilation, a standard clock with 
a measurable timescale is required. However, the challenge with objects such as 
quasars, and other cosmological sources such as gamma-ray bursts, is the com- 
plexity of the physical processes driving their variability. For quasars, where 
variability arises in the stochastic processes in the relativistic disk orbiting 
a supermassive black hole, the resultant luminosity fluctuations could poten- 
tially be dependent on a multitude of physical properties, including the mass 
of the central black hole, the degree of accretion, and the wavelength of the 
observations. 

To address this, the sample of quasars under consideration here was split 
into a number of subsamples of objects with similar intrinsic properties in 
terms of their bolometric luminosity and the rest wavelength of observations. 
The observations under consideration were taken in the (g,r,i) wavebands, 
and for the purpose of this study, the rest wavelength in each of the observed 
wavebands is defined to be 


X _ 4720A r  6415A y, — 7835A (1) 
a pee "O I+? * L+z2 


where z is the redshift of the quasar under consideration and the numerical 
values are representative of the observed wavebands. 

The quasar subsamples are presented graphically in Figure 2, which 
presents the rest wavelengths of each quasar in the (g,r,i) bands versus their 
bolometric luminosity. Each quasar has been colour-coded with its variability 
timescale, TpRw, assessed by fitting each observed light curve in each band 
with a damped random walk [see 27, for more detail]. Underlying the quasar 
sample are the regions of twelve subsamples under consideration in the colour 
salmon. These were chosen to have a width in rest wavelength and bolometric 
luminosity of AX = 1000A and A log(Lga/Lo) = 0.5. Note that the regions 
are continuous, with the top-left subsample spanning A\ = 900A > 1900A, 
and A log(Lgoi/Lo) = 46.7 — 47.2; the details of the subsamples are given in 
Table 2. From Figure 2 it is clear that these subsamples encompass the major- 
ity of quasars presented in this survey, and note that the combination of the 
three wavebands means that each subsample of quasars contains a broader 
distribution of redshifts than if the wavebands were considered individually. 
Hence this combination provides a redshift lever arm which constrains the 
presence of cosmological time dilation in each subsample. We also note that 
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a by-eye examination of Figure 2 is suggestive of a gradient in the variability 
timescale over the sample. 

Given that in each subsample the quasars possess similar rest wavelengths 
and bolometric luminosities, we make the assumption that they also possess the 
same characteristic intrinsic timescales, and hence any difference in timescale 
for a particular quasar subsample is due to the influence of cosmic time dilation 
and will show the appropriate dependence upon redshift. Of course, the physics 
of quasar variability is likely to depend on a number of factors, and so this 
assumption is considering that these quasars will exhibit similar variability 
properties in the mean; we discuss this point again at the conclusion of this 
study. 

For each quasar subsample (labelled with k), we model the observed 
variability timescales (i.e. logio(tp Rw /days)) as 


My, = Cy + nlogyg(1 + z) (2) 


where z is the redshift of the quasar and n is the power of the expected 
cosmological term, that is (1+ z)”. This model represents the variability with 
a different normalisation term, Ck, for each wavelength-luminosity bin, but 
demands the same cosmological dependence in terms of redshift. 

For the five distinct hypotheses considered, the normalisation terms, Ck, 
were allowed to vary, and so for the cases where n is considered to be a 
fixed value, this corresponds to an exploration of a twelve-dimensional pos- 
terior probability distribution. Physically, this situation reflects the situation 
where each wavelength-luminosity bin has a differing characteristic timescale, 
but a redshift dependence dependent upon the chosen value of n. For the 
remaining hypotheses, H2, where Ck and n are treated as free parameters, this 
corresponds to a thirteen-dimensional posterior probability distribution to be 
explored. 

To calculate the Bayesian evidence (also known as marginal likelihood) for 
each of these hypotheses, it is necessary to define a likelihood. It is impor- 
tant to note that the presented measurements and uncertainties of Tprw 
(in log;, space) are not symmetrical. Hence we represented the probability 
of each distribution of log, (tp rw /days) as a skewed Gaussian, specifically 
scipy.stats.skewnorm in the numerical approach which is written in python 
(represented as SN’). The 16%", 50° and 84!” percentiles of this skewed dis- 
tribution were fitted to the given values via a straight-forward optimization. 
This resulted in uncertainties of these percentiles of typically less than 2—3%. 
Hence we can define the log of the likelihood as 


loglL=S> X SSN logPDF(Mx, 0im(l,m € k)) (3) 


k l=g,r,i m=1...Nq 


where k sums over each of the subsample regions, l over the observed wave- 
bands and m over the number of quasars, Ng, in the sample. Also, 6), are the 
parameters for the skewed Gaussian representing the probability distribution 
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for log\9(tprw/days) for each of the quasars in each of the waveband, and 
l,m € k implies that the quasar should only be considered if its rest frame 
properties place it in subsample k. 

To explore the posterior probability distribution and calculate the Bayesian 
evidence (also known as marginal likelihood), we employed Diffusive Nested 
Sampling DNest4 [37], a variant of the nested sampling technique [38]. This 
allows for correct posterior sampling and marginal likelihood estimation even in 
the case where the constrained prior distributions are difficult to sample from 
or explore with Markov Chain Monte Carlo. For simplicity, we use uniform 
priors over the normalisation parameters, Ck, between 1 and 5, and for Ho 
where n is treated as a free parameter, a uniform distribution over n is adopted 
between -1 and 3; the posterior distribution of n for this hypothesis is presented 
in Figure 1. The normalisations, Ck, are well constrained and are reproduced 
for completeness in Figure 3 for H2. 


4 Conclusions 


This paper presents the detection of the cosmological dependence of the 
time dilation in a recent sample of almost two hundred quasars. These were 
monitored in multiple wavebands over a two-decade period, allowing the 
determination of a characteristic timescale by treating the observed quasar 
variability as a damped random walk. 

Through an assessment of the Bayesian evidence, it was found that the 
hypothesis considering the expected (1 + z) cosmological dependence provides 
a significantly better description of the data than the case where there is no 
dependence on redshift. In considering the redshift dependence of quasar vari- 
ability to be of the form (1 + z)”, where n is treated as a free parameter, the 
posterior distribution is found to be n = 1.28155, again consistent with the 
expected cosmological expansion of space. This detection of the cosmic expan- 
sion directly imprinted onto the variability of quasars further demonstrates 
that their observed properties are consistent with them being luminous and 
variable sources at cosmological distances, and counters previous claims that 
quasar variability is not intrinsic, but instead is due to external influences or 
non-standard physics. This has an immediate impact on various claims, such 
as the presence of a cosmologically significant population of microlensing black 
holes [e.g. 18, 22] or more esoteric ideas about the framework of the universe 
[39], and is further evidence that we inhabit an expanding relativistic universe. 

We do note that our result of n = 1.28%9'35 could be consistent with an 
offset from the expected cosmological value of n = 1 and could potentially 
indicate the presence of additional factors such as an evolution of quasars 
over cosmic time in addition to the time dilation due to cosmic expansion. 
Of course, we could imagine that quasar evolution over cosmic time could be 
responsible for the observed redshift dependence of the DRW time scale, but 
as we are considering similar quasars in terms of the bolometric luminosity and 
observed rest wavelength, it would be a curious coincidence for this evolution 
to result in a (1 + z) dependence to spoof cosmic expansion. Furthermore, if 
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quasar evolution were solely responsible for the observed DRW properties then 
the resulting lack of the expected cosmic time dilation would present a severe 
challenge to our cosmological model. However, it is important to note that 
there are some potential correlations of the DRW timescales with the inferred 
intrinsic properties of the quasars [e.g. 33], although these are not strong, and 
more extensive photometric datasets in terms of the number of quasars and the 
duration of their photometric lightcurves will be required to cleanly separate 
the influence of cosmic expansion from quasar evolution. 

In closing, we note that the lack of detection of the time dilation of quasar 
variability in previous studies is potentially due to the relatively small sample 
size in terms of the number of quasars under consideration [26], or the cadence 
of data sampling and characterisation of the quasar variability [19, 20]. Built 
on the observations of Stone et al. (2022) [27], this present study has demon- 
strated that we are now in an epoch where we have observations of a sufficiently 
large number of quasars spanning a broad range in redshifts, and observed over 
extended periods and with a cadence that overcomes their stochastic nature 
and results in an accurate characterisation of their variability, yielding a robust 
determination of the imprint of cosmological expansion on their light curves. 

Furthermore, with upcoming programs such as the Vera Rubin Observatory 

Legacy Survey of Space and Time (LSST), the number of quasars observed 
at high temporal cadence will rapidly increase and the measurement of cos- 
mological time dilation, and potentially the influence of quasar evolution, will 
become readily observable [e.g. 32]. 


Data Availability. The source data for this project is available at 
https://zenodo.org/record/ 5842449#.YipOg-jMJPY, with the details of 
the available FITS tables presented in Stone et al. (2022) [27]. Note that 
a revised version of this catalogue was recently released due to an error 
in some rest frame quantities. This revision does not impact any of the 
research presented in this paper. The software for this project is available at 
https://github.com/eggplantbren/QuasarTimeDilation. 


Code Availability. This project made use of several publicly available soft- 
ware packages, especially DNest4 [37] to undertake the exploration of the 
posterior probability space and calculate the Bayesian evidence by integrat- 
ing across this space. Further software packages employed include matplotlib 
[40], numpy [41], scipy [42]. Initial explorations of the posterior probability 
space were undertaken with emcee with corner plots prepared with corner 
[43]. The software employed as part of this project will be made available on 
reasonable request to the corresponding authors. 


Acknowledgements. We thank Stone et al. (2022) [27] for making their 
data and the results of their analysis publicly available. We also thank Scott 
Croom for his input and advice on quasar variability surveys. We further thank 
the teams responsible for creating and maintaining the various software pack- 
ages, detailed below, that this study has employed. GFL would like to thank 


Springer Nature 2021 TeX template 


8 On the Cosmological Time Dilation of Quasars 


the hospitality of the Lowell Observatory where the last stages of this work 
were completed during a period of isolation due to the contraction of covid. 


Author Contribution Statement. The project was conceived by GFL, 
including an initial exploration of the data, the definition of the models and 
hypotheses considered, the likelihood function and sampling of the posterior 
space. BJB undertook detailed sampling and calculating the Bayesian evidence 
using DNest4. Both authors discussed the results of the exploration in detail 
and determined the resulting conclusion. Both were responsible for the writing 
of the manuscript. 


Competing Interests Statement. The authors declare no competing 
interests. 


Tables. 


Springer Nature 2021 TeX template 


On the Cosmological Time Dilation of Quasars 


Table 1 The marginal likelihoods for the 
various hypotheses considered in this paper. 
As described in more detail in the Methods 
(Section 3), these were calculated with the 
diffusive nested sampling approach DNest4 [37]. 


Hypothesis n log Z 2 / Zmz 


Ho 0 -366.12 9.3 x 107-6 
Hy 1 -354.53 1 
H2 Free -356.52 0.14 
H3 -1 -390.13 3.5 x 10716 


Ha 2 -358.36 2.2 x 1072 
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Table 2 The properties of the survey subsamples presented in Figure 2, with the 
boundaries of the subsamples given by AX, rest wavelength, and A log; 9(LBoi/Lo), 
bolometric luminosity. The remaining columns give the number of quasar light curves in 
each subsample, Nqs, as well as the redshift range of those quasars, Az, and timescale for 
the observed variability as given by treating this as a damped random walk, 


A logio(to Rw /days). 

Subsample A (A) Alog(Leoi/Lo) Nas Az A log19(tp Rw /days) 
1 900 — 1900 46.7 — 47.2 37 1.60 — 4.15 2.68 — 4.03 
2 1900 — 2900 46.7 — 47.2 27 0.81 — 3.00 2.83 — 4.11 
3 900 — 1900 46.2 — 46.7 74 1.55 + 3.98 2.61 — 4.23 
4 1900 — 2900 46.2 — 46.7 111 1.11 > 3.01 2.49 — 4.42 
5 2900 — 3900 46.2 — 46.7 22 1.11 — 1.70 3.03 —> 3.93 
6 900 — 1900 45.7 — 46.2 30 1.48 — 2.80 2.70 —> 3.71 
7 1900 — 2900 45.7 — 46.2 101 0.68 — 2.80 2.55 — 3.92 
8 2900 — 3900 45.7 — 46.2 58 0.60 — 1.69 2.63 — 4.03 
9 3900 — 4900 45.7 — 46.2 11 0.60 — 1.00 2.66 — 3.84 
10 1900 — 2900 45.2 — 46.7 27 0.63 —> 1.45 2.30 — 4.30 
11 2900 — 3900 45.2 — 46.7 31 0.47 — 1.45 2.27 — 4.30 
12 3900 — 4900 45.2 — 46.7 20 0.47 — 0.98 2.33 — 4.20 
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Figures. 


W126) 55 


—1.0 —0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 
n 


Fig. 1 The posterior distribution of n, where the redshift dependence of the observed time 
dilation is given by (1 + z)”, for the Bayesian exploration of H2, where this index is treated 


as a free parameter in the analysis. From this distribution, n = 1.2870 38, where the best 


fit value is taken as the median (50t percentile), whilst the uncertainties represent the 
16t? and 84t} percentiles. This distribution was determined through an exploration of the 
posterior probability space with DNest4 [37]. 
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Fig. 2 The entire quasar sample under consideration as a function of rest wavelength and 
bolometric luminosity, colour-coded with the DRW timescale, Tp>Rw. The underlying rect- 
angles in salmon pink represent the boundaries of the subsamples employed in the analysis 
presented in this paper. Inset is the labelled numbers of the fields. 


Springer Nature 2021 {TX template 


On the Cosmological Time Dilation of Quasars 13 


46.7 > 47.2 


45.7 > 46.2 


2.0 2.5 3.0 3.5 2.0 2.5 3.0 3.5 
Ck 


Fig. 3 The posterior distributions for the normalisation parameters, Ck, for #2 where n is 
treated as a free parameter. The bolometric luminosity range for each of the normalisation 
parameters is given in the upper left of each panel (see Table 2). As with Figure 1, these 
were the result of the sampling of the posterior probability distribution with DNest4. 
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